
## FIGURES 3 AND 4

# load, combine data
library(rio)
library(tidyverse)

# estimate models
library(brglm2)

# cluster standard errors
library(lmtest)
library(sandwich)

# generate plots
library(sjPlot)
library(scales)
library(ggpubr)

# quantify effects
library(ggeffects)




# LOAD DATA
glm_df <- import("glm_data.csv") %>%
  filter(iso3c!="TTO" & iso3c!="GUY" & iso3c!="SUR") %>%
  filter(!is.na(approval) & !is.na(opp1vote) & !is.na(number_protests))



#### PREDICTED PROBABILITIES

# based on Model 6, Table 2
brglm6_plot <- glm(doc ~ approval + opp1vote + number_protests + previous_doc + duck + partyage + polity2 + left_exec + any_election +
                     resource_rents_lag + discovery + real_oil_price_log + gdp_pc_constant_lag + gdp_growth_lag + imf_program + crisis +
                     #as.factor(iso3c) + 
                     quarter + year, 
                   family = binomial(link="logit"), 
                   method = "brglmFit",
                   data = glm_df)


plot1 <- plot_model(brglm6_plot, type = "pred", terms = c("approval[all]"),
                    se = vcovHC(brglm6_plot, type = 'HC0', cluster = c('iso3c', 'Summary'))) +
  theme_minimal() + labs(x = "Executive Approval", y = expression(paste("Predicted Probability of ", italic("Any Document"))),
                         title = " ") + scale_x_continuous(labels = label_percent(scale = 1)) + 
  scale_y_continuous(lim = c(0,0.3), labels = label_percent(scale = 100, accuracy = 1))
plot2 <- plot_model(brglm6_plot, type = "pred", terms = "opp1vote[all]",
                    se = vcovHC(brglm6_plot, type = 'HC0', cluster = c('iso3c', 'Summary'))) +
  theme_minimal() + labs(x = "Opposition Vote Share", y = expression(paste("Predicted Probability of ", italic("Any Document"))),
                         title = " ") + scale_x_continuous(labels = label_percent(scale = 1)) + 
  scale_y_continuous(lim = c(0,0.15), labels = label_percent(scale = 100, accuracy = 1))
plot3 <- plot_model(brglm6_plot, type = "pred", terms = "number_protests[all]",
                    se = vcovHC(brglm6_plot, type = 'HC0', cluster = c('iso3c', 'Summary'))) +
  theme_minimal() + labs(x = "Number of Protests", y = expression(paste("Predicted Probability of ", italic("Any Dcument"))),
                         title = " ") + 
  scale_y_continuous(lim = c(0,0.9), labels = label_percent(scale = 100, accuracy = 1))


brglm12_plot <- glm(fund_doc ~ approval + opp1vote + number_protests + previous_fund_doc + duck + partyage + polity2 + left_exec + any_election +
                      resource_rents_lag + discovery + real_oil_price_log + gdp_pc_constant_lag + gdp_growth_lag + imf_program + crisis +
                      #as.factor(iso3c) + 
                      quarter + year, 
                    family = binomial(link="logit"), 
                    method = "brglmFit",
                    data = glm_df)

plot4 <- plot_model(brglm12_plot, type = "pred", terms = "approval[all]",
                    se = vcovHC(brglm12_plot, type = 'HC0', cluster = c('iso3c', 'Summary'))) +
  theme_minimal() + labs(x = "Executive Approval", y = expression(paste("Predicted Probability of ", italic("Fund Document"))),
                         title = " ") + scale_x_continuous(labels = label_percent(scale = 1)) + 
  scale_y_continuous(lim = c(0,0.3), labels = label_percent(scale = 100, accuracy = 1))
plot5 <- plot_model(brglm12_plot, type = "pred", terms = "opp1vote[all]",
                    se = vcovHC(brglm12_plot, type = 'HC0', cluster = c('iso3c', 'Summary'))) +
  theme_minimal() + labs(x = "Opposition Vote Share", y = expression(paste("Predicted Probability of ", italic("Fund Document"))),
                         title = " ") + scale_x_continuous(labels = label_percent(scale = 1)) + 
  scale_y_continuous(lim = c(0,0.15), labels = label_percent(scale = 100, accuracy = 1))
plot6 <- plot_model(brglm12_plot, type = "pred", terms = "number_protests[all]",
                    se = vcovHC(brglm12_plot, type = 'HC0', cluster = c('iso3c', 'Summary'))) +
  theme_minimal() + labs(x = "Number of Protests", y = expression(paste("Predicted Probability of ", italic("Fund Document"))),
                         title = " ") + 
  scale_y_continuous(lim = c(0,0.9), labels = label_percent(scale = 100, accuracy = 1))


# ggsave("figure3.pdf", ggarrange(plot1, plot4), width = 9, height = 4)
# ggsave("figure4-1.pdf", ggarrange(plot2, plot5), width = 9, height = 4)
# ggsave("figure4-2.pdf", ggarrange(plot3, plot6), width = 9, height = 4)


# some concrete numbers to include in the results discussion

# effect of executive approval
ggpredict(brglm6_plot, terms = "approval[82.45]")
ggpredict(brglm12_plot, terms = "approval[82.45]")

# effect of opposition vote share
ggpredict(brglm6_plot, terms = "opp1vote[52.2]")
ggpredict(brglm12_plot, terms = "opp1vote[52.2]")

# effect of number of protests
ggpredict(brglm6_plot, terms = "number_protests")
ggpredict(brglm12_plot, terms = "number_protests")
